An algorithmic approach to simulate Hamiltonian dynamics 
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We propose an iterative algorithm to simulate any n-qubit Hamiltonian dynamics. The simulation 
entails decomposing a unitary operation into a product of simpler unitary operations, where each 
unitary operation can be easily implemented in a quantum simulator of choice. Our algorithm 
exploits the completeness of the chosen operator basis to intuitively understand the structure of the 
system dynamics, and scaling properties of the simulation. We illustrate the algorithm by simulating 
a unitary that transports quantum information from one end of a n-qubit chain to the other driven 
by an Ising Hamiltonian (an operation known as quantum state transfer (QST)). This is applied to 
provide the first experimental verification of the QST protocol by Di Franco et. al. (Phys. Rev. 
Lett. 101, 230502 (2008)) in an NMR quantum simulator. We also experimentally test the resilience 
of this protocol to static errors in the coupling parameters of the simulated Hamiltonian. 

PACS numbers: 03.67.Ac, 76.60.-k, 03.67.Lx 



Introduction: Feynman [l| alluded to the fact that one 
could exploit quantum mechanical evolution to simulate 
the dynamics of another quantum system exponentially 
faster than in any classical computer. This has been one 
of the key motivations to build a quantum information 
processor (QIP) 0). Simulation not only allows one to 
study the dynamics of another system, perhaps not nat- 
urally available, but also to execute a certain computa- 
tional task (an algorithm) Q. More precisely, the sim- 
ulation problem is the desire to efficiently implement an 
unitary operation U — exp(-i'Ht) on ones choice of QIP, 
where H is an n-qubit Hamiltonian. An efficient simu- 
lation is one which requires physical resources that scale 
polynomially with the problem size n. 

The Hamiltonian of any QIP consists of the sum of the 
intrinsic Hamiltonian "Hint , and a part which depends on 
parameters controlled experimentally by externally ap- 
plied fields Hext(t) Their specific form, of course, 
depends upon the inherent physical nature of the QIP. 
The parameters are chosen such that the following holds 
to an arbitrary precision: 



U « U sim = Texp 



dt{n mt + n cxt (t)) 



(1) 



where t is a finite time interval. If a QIP allows universal 
control, i.c in principle it can simulate any desired uni- 
tary operation, then there exist the control parameters 
such that Eq.m holds @|. However, finding them is the 
real challenge [6]. Instead, t/ s i m is obtained as a sequence 
(product) of unitaries that are easily implementable in 
the QIP: U s - lm — XJ\Ui....U m [7]. In general m scales ex- 
ponentially with the number of qubits n, unless U has 
additional structure, for instance, its matrix representa- 
tion is sparse, or it has certain kinds of symmetry with 
respect to its action on the space of qubits |8|] . Hence, the 
product decomposition of U unravels the physical struc- 
ture of U: how it entangles different qubits, which part 



of the Hilbert space it accesses and how the resources re- 
quired for its simulation scale with the number of qubits. 

Unfortunately, the conventional method to obtain the 
decomposition does not allow one to probe such struc- 
ture of U . In the standard methods @ (for small n) the 
decomposition is achieved by resolving r into finer time 
steps Stj, with the assumption that Stj — > 0. A first 
order approximation is invoked, and numerical optimiza- 
tion techniques are employed: A guess decomposition U g 
is made, which is then compared with the desired U via 
an overlap function (like fidelity); U g is then bought ar- 
bitrarily close to U by iteration [9(. Since the process 
requires infinitesimal time steps Stj, the decomposition 
obtained is discontinuous and highly nonintuitive. In this 
letter, we rectify this problem by algorithmically decom- 
posing U instead into finite time step buil ding blocks that 
are elements of a complete operator basis [lfj . Moreover, 
the decomposition obtained by our method for small n 
may allow one to deduce the same for arbitrary n. 

In this letter we demonstrate the utility of our algo- 
rithm by explicitly determining the decomposition of the 
unitary that causes quantum state transfer (QST) in a 
linear n-spin chain [111 ]. QST physically refers to the 
transport of quantum information from one end of the 
chain to the other; allowing spin chains to act as ana- 
logues of wires in spin-based QIP architectures [l2| ■ QST 
is slated to play an important role in the initialization of 
quantum computers, and in linking together the quan- 
tum registers of two different spin based QIPs ■ While 
QST can be achieved by a hopping transfer between suc- 
cessive spins of the chain, this requires locally addressing 
each qubit, and is not experimentally viable for large n. 
Instead, a global control field on the spin chain, which 
is in general entangling, maybe more efficient [111 ]. For 
example, QST can occur if the spins are driven by an 
Ising Hamiltonian T-Li (Eq. [3] below). But this requires 
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manufacturing a chain of spins coupled to each other in 
a very specific way - a task that is fraught with enor- 
mous difficulty 13]. So, while several protocols have been 
suggested for QST using global control Hamiltonians like 
Hi, there is no experimental means to test their efficiency 
or their robustness to fabrication errors in the chain. Ex- 
perimentally simulating the quantum chain allows one to 
probe such issues. 

Here, we apply the decomposition to experimentally 
simulate QST in an Nuclear Magnetic Resonance (NMR) 
QIP [14j . We implement the QST protocol recently pro- 
posed by Di Franco et al [l5[ that employs Hi, and does 
not require the spin chain to be initialized. We find that 
the resources for simulation scale linearly with the num- 
ber of qubits, and that the Ising Hamiltonian entangles 
all the qubits of the chain. The algorithm is also used 
to obtain decompositions to experimentally test the re- 
silience of this QST protocol to errors in the static cou- 
plings between spins. 

The Quantum State Transfer. Given a n-qubit chain, 
initially in the statel Vi) ' ' ' iV'n)) Di Franco et al sug- 
gested a protocol [15[ by which QST can be achieved 
by evolving the spin chain under the nearest-neighbour 
Ising Hamiltonian Hi (given below) for a time to — ir/AJ 
■ Uqst = exp(—it Hi). The importance of this protocol 
is that it bypasses the initialization of the spin medium 
(except the last qubit); hence the core of the spin chain 
is not under local control [15j . 

Uqsr\1>i)\'h---il> n -i)\0) (2) 

= -L [|o>|v„-i ■ -ihMi) + Wlri-i • ••W-)(ixhk»] 

where (V'i'h/'i) = OVi A state locally equivalent to \ipi) 
can be recovered at the end of the chain by measuring 
the first qubit in the Z basis. 



The Ising Hamiltonian Hi has the form [15j . 



Hi = 22 JjI 3Z I (j+1)z + ^2 B i I 'j^ 

3=1 3=1 



(3) 



J's are a spin-1/2 Pauli matrices and we assume stan- 
dard commutation relations hold. J,- refers to the inter- 



neighbour coupling strength, Jj = 2J^Aj(n — j); while 
Bj refers to the coupling strength of a spin site with a 
local magnetic field, Bj = 2 J y / {2j - l)(2n - 2j + 1). 

The algorithm: In this letter we illustrate the algorithm 
by explicitly obtaining a decomposition of the 3-qubit 
QST unitary: Uqst = exp (— itoHi), for to — 7r/4J. Fol- 
lowing the steps outlined here, the algorithm can be eas- 
ily generalized for any unitary U. We will decompose 
Uqst in the well-known NMR product operator basis 
.17]. 



where {a, /3, 7} € {x,y,z}. The aim of the algorithm is 
to find the m basis-elements Dk(& B) and angles 9 k such 
that, 



f^QST = Yi exp(-i0 k D k ) , 



(5) 



fc=i 



Since the D k s act on different parts of the qubit space, 
the decomposition reflects the internal structure of Uqst 
in the chosen basis. 

Determining D k s in Eq. ([5]) can be looked upon as 
a search problem, where the search-space is the set B. 
For an n-qubit system this space is of exponential size 
2 2n . The central step in the algorithm is to obtain the 
decomposition of Uqst by iteratively reducing this search 
space (as explained below). By using exp(— i^-Dfc) = 
(cos -y 1 — i sin %-D/.) , Eq. ([5]) can be re-expressed as, 

Uqst = C 1 l + C Dl D 1 +"- + C Dm D rn + C DlIh D 1 D 2 + 



I, 



C Dl ... Dm D 1 ---D m , (6) 



which can be simplified by using the matrix product oper- 
ation as the expansion of the unitary in the complete ba- 
sis B. In the special case that Go = {1, D\,D 2 , • • • , D m } 
forms a group, Eq. © reduces to, 



Uqst = Ci1+Cd 1 Di- 



Cn D ri 



fe=i 

(7) 

where the coefficients C Dk = Tr(U QST Dl)/Tr(DlD k ). 
Note that Go can always be formed by suitable inclusion 
of elements from B. The closure of the group Go implies 
that in this case all the elements D k in the decomposition 
ofEq. (O belong to Go . To see this, consider for example 
the explicit expansion of the Uqst as a sum of elements 
in B: 

C^qst = —= [1 — ihx + 2Ilxhx + %hyhy + 2Ilzhz 



2V2 



i±hzh*hz] (8) 



lAI\ v l2 x I 3 y 

The elements that figure in this basis expansion form a 
group, 

Go = {l,hx,2h x^3x 1 x *2x ^3x 1 



^lylzxHy 7 4-/i 2 /2;E-f3z}, 



(9) 



Hence, Uqst can be decomposed (as in Eq. ([5])) by using 
D k s that exclusively belong to Go, i.e. the search space 
for finding D k is now Go C B. The closure of the group 
also provides a means to iteratively reduce this search 
space: by extinguishing part of Uqst to get J7qg T so that 
the only elements in its basis-expansion (similar to Eq. 
([7])) now form a subgroup of Gq. For instance, consider 



the unitary operator Uqq T , 



U, 



(i) 

QST 



B — {l,Ij a ,2Ij a Ikp,4:Ij a IkaIlj} 



(4) 



Uqst x exp (-i^<iI ly I 2x I 3y J 
hi- il 2x + Zhzhz + i^hzhxhz 



(10) 
(11) 



3 



Now, since G% = {l.hx^hzhzAhzhxhz} is a sub- 
group of Go; the search-space is reduced to Gi c Go C B. 
A similar iterative reduction can be continued by extin- 
guishing part of Uqsi to get C^Q 2 g T , so that the search- 
space G2 is a subgroup of Gi , and so on. More generally, 
let us assume that forming Uqq t from Uqst requires ml 
elements, 



C/qst x JJ cxp(i6 k D k 



= U, 



(i) 

QST- 



(12) 



. k=l 



These m! basis-elements are precisely those that figure 
in the decomposition of E/qst in Eq. ([5]). We define the 
norm of the space spanned by the elements of a set G 
(G C B) in a unitary U as, 



M G {U) 



k 



C D J 2 where, D k e G 



With this definition, 

M?i(Cqst) + A/( Go _ Gi )(C/qst) 



(13) 



(14) 



The role of the ml elements D k in Eq. (fl"2|) is to make 

■^(G -Gi)(^qst) = 0- They way they are ordered (in the 
decomposition) is chosen by dynamic programming (ljj : 
at each of ml steps one affects a maximum reduction in 



norm; so that at the end of m'-steps, the norm of the 
space spanned by (Go — Gi) in t/qg T is arbitrarily close 
to 0. For this process, it is sufficient that the D k 's belong 
exclusively to (Go — Gi). The angles 9 k in Eq. (TT2")) 
are chosen in such a way that the associated D k has no 
support in Uqg T . Recurrence relations for the dynamic 
programming process and their scaling will be presented 
elsewhere [161 ]. In summary, the 3-qubit Uqst has the 
decomposition: 

Uqst = exp {i^hyhxh^j exp (i^hzhxhz^ exp 

(15) 

The advantage of our algorithmic approach is that by ob- 
taining the decomposition of a unitary for small n, one 
maybe able to deduce its decomposition for arbitrary n. 
For J7qst, this decomposition is given in Eq. (|16[) . where 
the upper equation is for odd n, and the bottom one for 
even n. The number of basis elements in the decomposi- 
tion (hence the resources required for simulating C/qst) 
scale linearly with n. Since, in general, directed informa- 
tion transport along a quantum chain would involve all 
the qubits, at least 0(n) resources are needed to simulate 
it. Our decomposition is hence minimal in this regard. 
Eq. (fi"6|) also reflects the fact that C/qst entangles all 
spins in the chain as suggested by the protocol in Eq. 
©• 



L«/2J 



C/qst = JJ ex P ( fy n ~ 2k+ll k y I{n-k+i) v (g) I 3 * 

j=k+i 



exp i—2 



fc=i 



1 — k 



n—2k+l t t 

J-kz-l(n-k+l)z 



7 1x1 cx P(" z f J ("+i)/ 2 



j=k+l 



[161 



Experimental Simulation of Ising Hamiltonian: NMR 
is a pioneering architecture for experimental QIP, where 
the development of sophisticated control methods over 
several decades has led to the exploitation of the full 
range of the system dynamics 0] . In an NMR QIP 17 1 , 



Htm = I E"=i Ijlj* + I E™=2 J2k< 3 JkjhJjz, where 7i 
is the gyromagnetic ratio of the jth qubit (spin), and 
J k j is the J-coupling between spins j and k. The exter- 
nal control field H e xt(t) = u x(t) J2j Ijx + u y (t) J2 3 Ijy is 
implemented via radio frequency (RF) pulses (the pulse 
sequence) that generate single qubit rotations about the 
x/y-eods — this ensures universal control 

In our experiments, we use 13 C labeled 13 CHFBr2, 
where 1 H, 13 C and 19 F are the three-qubits. Experiments 
were performed at room temperature in a 11. 7T magnetic 
field on an AV-500 spectrometer with a triple resonance 
QXI probe. The 1 H, 13 C and 19 F resonance frequen- 
cies at this field are 500MHz, 470MHz and 125MHz re- 
spectively. The couplings between the three-qubits are 



Jhc = 224.5Hz, J HF = 49.7Hz and J FC = -310.9Hz. 

Usin g hard p ulses alone or by standard refocusing tech- 
niques [17J, ll9j , each element of the product operator ba- 
sis B (Eq. (j3J) can be implemented experimentally. One 
can arrive at the pulse sequence that implements [/qst 
by piecing together the pulse sequences to construct all 
the elements in the decomposition. Moreover, in Eq. f| 15[) 
(and Eq. (1TB)) ) all elements in the decomposition can be 
implemented by only using nearest neighbour couplings. 
The other couplings are refocused during the net period 
of the sequence. Hence an effective qubit chain can be 
created in the NMR QIP, and the pulse sequence alters 
the total NMR Hamiltonian so that it evolution is iden- 
tical to that of the Ising Hamiltonian Hi after tg. We 
will employ 1 H, 13 C and 19 F as a 3-qubit chain in that 
order; the Jhf coupling is refocused and not used. 

Let \ipi) = a\0) + j3\l) be the state of the first qubit 
which is to be transferred along the spin chain. The 
protocol 15[ requires that the last qubit be initialized 
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Sequence for U QST —cxp(-iw'H I /4J) 



1 H 



1.57 t 1.57 
2nJ l2 | 2nJ 2:t 



1.57 

27T./..> 



1.57 
27rJ 12 



1 I 2»? K [ 

X . ~TJ X 



1.57 f 1.57 
2*jj3 1 27r,7 1; I 



1.57 



(a) 



(i) Verification of QST protocol (ii) V., with imperfect couplings 

l.Ofr- 



■a °- 5 - 

c 

0) - 

4-J 7 

£ o.o» 

a 

Oi -0.5- 

-1.0s- 



T 



— Theory 

• • Experiment 



tt/ 2 3tt/2 

Angle <p (rad) 




tt/2 ff 3tt/2 2tt 

Angle ^ (rad) 





=-1.756J 


• • i.7; =«, =0 


&7, =<S,7 2 =1.270,7 


♦ ♦ «j =W, 


= -0.756,7 


i — ^ SJ, =<5,7 2 =0.6677 





(b) 

Sequence for H r with imperfect couplings of SJ ] —5J 2 = —0.756.7 



, Q.14 
fl 2jtJ m 



1.57 
■ 2irJ M 



1.57 , 

! 27T./ la | 



13 C 



5.95 f 6.15 * 1.57 f 4.CJ7 f 1.57 f * 1.57 4.51 X 1.57 , 



27rJ J3 i i In J 12 i I W a |||; 2tJ 1b | 2^J Ba I 27TJ, 
X X X XX 

n irsti n n 



(c) 

Sequence for with imperfect couplings of fl,/, =0.667.7 



1 



3.46 3.46 
! 2ttJ 2 , ! ! 2tt,7 12 



1.57 -? 1.61) ^ 1.57 - ? 1.57 , 
; 2irJ u | 271-./,;, |j 2ttJ u [|; 2^ | 



27rJ ia [j; 2^ I 27rJ. : , I 2tvJ v , 



XX XX X X X X XX XX, on XXX 



, 4.89 ' 
2mZ» 



(d) 



FIG. 1: (Color online) (a) Pulse sequence for Uqst- Filled 
wide and narrow bars are 7r and 7r/2 pulses respectively. Un- 
filled bars represent pulses with other flip angles. Dotted lines 
are erstwhile pulses that have been removed by simplification 
and can be ignored. The numbers between pulses are de- 
lays. The black arrows show the sequences for the individual 
operators (boxes) in the decomposition of Uqst- (b)(i) Theo- 
retical (line) and experimental (points) normalized signal in- 
tensity on the third qubit, obtained as the average intensity of 
the spectral lines corresponding to the state ( 1 00) +i|ll))|V'/) 
(shown in the inset for if — n/2). Here fo = — 40.26kHz(-85.6 
ppm). (ii) Experimental signal intensity on the third qubit for 
different static errors in the Ising Hamiltonian couplings. Two 
representative pulse sequences used in (ii) are shown in (c), 
(d). The unfilled bars in (c) and (d) are 76.52° and 168.71° 
X pulses respectively. 



to |0), while the second qubit can be in any arbitrary 
state. For ease of measurement, we will fix the initial 
state of the second qubit also to |0). The initial state of 
the spin chain, |?/>/00), was prepared from a 1 000) psue- 
dopure state using a (p-angle X pulse on the first qubit, 
where <p = tan -1 — . The pulse sequence for Uqst (Fig 
1(a) ) was then applied. However, for the readout of 
projective measurement is not possible in an ensemble 
QIP architecture like NMR. Instead, to verify that QST 
has indeed occurred, we applied a CNOT(l,3) gate, forc- 
ing the system in the state ( 1 00) + £|ll))|^r). The state 
of the third qubit was measured with a Y pulse with the 
receiver coil in the X direction. The obtained NMR sig- 
nal shows excellent correlation with sin</j (Fig. l(b)[ ), 
showing that the state \ipj) has been transferred to the 
last qubit. 

How does the effectiveness of the QST protocol change 
if the couplings Jj of the Ising Hamiltonian are imper- 
fect? We consider the case where the deviations of the 
couplings from their Ising values are equal, i.e SJi = <5 J 2 ; 
and use the NMR simulator to test the QST process by 
determining sequences that simulate the unitaries result- 
ing from the changed Hamiltonians. Two representative 



examples of these sequences are shown in Fig. 1 (c) and 
1(d) Not surprisingly, from a decomposition point of 



view, these sequences have a larger structure compared 
to the Ising Hamiltonian with perfect couplings. Our 
experimental results (Fig. l(b)[ i) show that the state re- 
covered is corrupted by an additional phase, and the de- 
gree of corruption increases with the error of couplings. 
Consequently, the final state on the last qubit is not a 
faithful reflection of \tpi), and the QST process is imper- 
fect. Moreover, the transfer efficiency - characterized by 
the maximum signal amplitude - decreases with the error 
of couplings. 

Summarizing, we have presented an algorithm to pro- 
vide a decomposition for any desired unitary operator in 
a manner that provides an intuitive understanding of its 
internal structure, and which allows the unitary to be 
simulated by an NMR QIP. We used it to experimentally 
verify the quantum state transfer protocol [lol ] in spin 
chains based on the Ising Hamiltonian. 

This letter is dedicated to the memory of Jharana Rani 
Samal. 
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